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Abstract 

Smoothed Particle Hydrodynamics is a multidimensional Lagrangian method of nu- 
merical hydrodynamics that has been used to tackle a wide variety of problems 
in astrophysics. Here we develop the basic equations of the SPH scheme, and we 
discuss some of its numerical properties and limitations. As an illustration of typi- 
cal astrophysical applications, we discuss recent calculations of stellar interactions, 
including collisions between main sequence stars and the coalescence of compact 
binaries. 



1 Smoothed Particle Hydrodynamics 

Smoothed Particle Hydrodynamics (SPH) is a Lagrangian method that was 
introduced specifically to simulate self-gravitating fluids moving freely in three 
dimensions. The key idea of SPH is to calculate pressure gradient forces by 
kernel estimation, directly from the particle positions, rather than by finite 
differencing on a grid, as in older particle methods such as PIC. SPH was 
first introduced by Lucy (1977) and Gingold & Monaghan (1977), who used 
it to study dynamical fission instabilities in rapidly rotating stars. Since then, 
a wide variety of astrophysical fluid dynamics problems have been tackled 
using SPH (see Monaghan 1992 for an overview). In addition to the stellar 
interaction problems described in §2, these have included planet and star for- 
mation (Nelson et al. 1998; Burkert et al. 1997), supernova explosions (Herant 
h Benz 1992; Garcia-Senz et al. 1998), large-scale cosmological structure for- 
mation (Katz et al. 1996; Shapiro et al. 1996), and galaxy formation (Katz 
1992; Steinmetz 1996). 
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1.1 SPH from a Variational Principle 



A straightforward derivation of the basic SPH equations can be obtained from 
a Lagrangian formulation of hydrodynamics (Gingold & Monaghan 1982). 
Consider for simplicity the adiabatic evolution of an ideal fluid with equation 
of state 

P = Ap\ (1) 

where p is the pressure, p is the density, 7 is the adiabatic exponent, and 
A (assumed here to be constant in space and time) is related to the specific 
entropy (s oc In A). The Euler equations of motion, 

_ = _ + (B .V). = --V R (2) 

can be derived from a variational principle with the Lagrangian 

L = j{\v 2 -u[p{r)\}pd'x. (3) 



Here u[p] = p/[(j — l)p] = Ap 1 1 /(j — 1) is the specific internal energy of the 
fluid. 

The basic idea in SPH is to use the discrete representation 



N 

L S PH = J2 m i 
i=l 



(4) 



for the Lagrangian, where the sum is over a large but discrete number of small 
fluid elements, or "particles," covering the volume of the fluid. Here mi is the 
mass and Vi is the velocity of the particle with position fj. For expression (4) 
to become the Lagrangian of a system with a finite number N of degrees of 
freedom, we need a prescription to compute the density pi at the position of 
any given particle i, as a function of the masses and positions of neighboring 
particles. 

In SPH, the density at any position is typically calculated as the local average 

p(f) = ^2m j W(f-f j ;h), (5) 
3 



where W(x; h) is an interpolation, or smoothing, kernel of width ~ h. Nec- 
essary constraints on the kernel W(x; h) are that (i) it integrates to unity 
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(consequently the integral of eq. (5) over all space automatically gives the 
total mass of the system), and (ii) it approaches the Dirac delta function S(x) 
in the limit h — > 0. 

Equation (5) gives, in particular, the density in the vicinity of particle i as 
Pi = pifi), and we can now obtain the equations of motion for all the particles. 
Deriving the Euler-Lagrange equations from Lsph we get 




where Wij = W(fi — fj] h) and we have assumed that the form of W is such 
that = Wji. The expression on the right-hand side of eq. (6) is a sum 
over neighboring particles (within a distance ~ h of fi) representing a discrete 
approximation to the pressure gradient force [— (l/p)Vp]j acting on particle i. 

The following energy and momentum conservation laws are satisfied exactly 
by the simple SPH equations of motion given above 

and 




where Ui = Pi/ [('J — l)pi] . Note that energy and momentum conservation in 
this simple version of SPH is independent of the number of particles N. 

Typically, a full implementation of SPH for astrophysical problems will add 
to eq. (6) a treatment of self-gravity (e.g., using one of the many grid-based 
or tree-based algorithms developed for N-body simulations) and an artificial 
viscosity term to allow for entropy production in shocks. In addition, we have 
assumed here that the smoothing length h is constant in time and the same for 
all particles. In practice, individual and time- varying smoothing lengths hi(t) 
are almost always used, so that the local spatial resolution can be adapted to 
the (time-varying) density of SPH particles (see Nelson & Papaloizou 1994 for 
a rigorous derivation of the equations of motion from a variational principle 
in this case). Other derivations of the SPH equations, based on the applica- 
tion of smoothing operators to the fluid equations (and without the use of a 
variational principle), are also possible (see, e.g., Hernquist & Katz 1989). 
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1.2 Basic SPH Equations 

In this section, we summarize the basic equations for various forms of the 
SPH scheme currently in use, incorporating gravity, artificial viscosity, and 
individual smoothing lengths. 

1.2.1 Density and Pressure 

The SPH estimate of the fluid density at r*j is calculated as p« = J2j m j^ij [ c f- 



eq. (5)]. Many recent implementations of SPH use a form for Wij proposed by 
Hernquist & Katz (1989), 



This choice guarantees symmetric weights = even between particles i 
and j with different smoothing lengths. For the interpolation kernel W(r; h), 
the cubic spline 



(Monaghan & Lattanzio 1985) is a common choice. Eq. (10) is sometimes 
called a "second-order accurate" kernel. Indeed, when the true density p(r) 
of the fluid is represented by an appropriate distribution of particle positions, 
masses, and smoothing lengths, one can show that pi = p(fi) + 0(hj) (see, 
e.g., Monaghan 1985). 

Depending on which thermodynamic evolution equation is integrated [see 
eqs. (26) and (27) below], particle % also carries either the parameter u h the 
internal energy per unit mass in the fluid at fj, or A^, the entropic variable, 
a function of the specific entropy in the fluid at fj. Although arbitrary equa- 
tions of state can be implemented in SPH, here, for simplicity, we consider 
only polytropic equations of state. The pressure pi at fi is therefore related to 
the density by 



^[Wdn-fjl-h^ + wdn-fjl-hj)}. 



(9) 




h — 



(10) 



Pi = (7 - 1) Pi u 



-1 1 



(11) 



or 



Pi = AipJ. 



(12) 
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The speed of sound in the fluid at f; is q = (7^/ 'Pi) 1 ^ 2 ■ 



1.2.2 Dynamical Equations and Gravity 
Particle positions are updated either by 

dfi 



a* *• (13) 



or the more general XSPH method 



— = Vi + e ^j— n Wn (14) 

m j Pij 



where pij = (pi + Pj)/2 and e is a constant parameter in the range < 
e < 1 (Monaghan 1989). Eq. (14), in contrast to eq. (13), changes particle 
positions at a rate closer to the local smoothed velocity. The XSPH method 
was originally proposed as a way to minimize spurious interparticle penetration 
across the interface of two colliding fluid streams. 

Generalizing equation (6) to account for gravitational forces and artificial vis- 
cosity (hereafter AV), the velocity of particle % is updated according to 

dVi _ _(Grav) -*(SPH) 



dt 



where aj Grm '' ) is the gravitational acceleration and 



-(SPH) 



= -E 



m, 



£ + £ +n- 

A + pV + 1 



(16) 



Various forms for the AV term are discussed below. The AV ensures that 
correct jump conditions are satisfied across (smoothed) shock fronts, while the 
rest of equation (16) represents one of many possible SPH-estimators for the 
acceleration due to the local pressure gradient (see, e.g., Monaghan 1985). 

To provide reasonable accuracy, an SPH code must solve the equations of 
motion of a large number of particles (typically N >> 1000). This rules out a 
direct summation method for calculating the gravitational field of the system, 
unless special purpose hardware such as the GRAPE is used (Steinmetz 1996; 
Klessen 1997). In most implementations of SPH, particle-mesh algorithms 
(Evrard 1988; Rasio & Shapiro 1992; Couchman et al. 1995) or tree-based 
algorithms (Hernquist & Katz 1989; Dave et al. 1997) are used to calculate 
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the gravitational accelerations a\ r . Tree-based algorithms perform better 
for problems involving large dynamic ranges in density, such as star formation 
and large-scale cosmological simulations. For stellar interaction problems like 
those described in §2, density contrasts rarely exceed a factor ~ 10 2 — 10 3 and 
in those cases grid-based algorithms and direct solvers are generally faster. 
Tree-based and grid-based algorithms are also used to calculate lists of nearest 
neighbors for each particle exactly as in gravitational iV-body simulations. 



1.2.3 Artificial Viscosity 

For the AV, a symmetrized version of the form proposed by Monaghan (1989) 
is often adopted, 

Pij 



where a and j3 are constant parameters, c^- = (q + Cj)/2, and 



(vi-Vj)-(n-rj) 



if (Hi - Vj) ■ (fi - fj) < 
if (vi - vj) ■ (f - fj) > 



(18) 



with hij = (hi + hj)/2. This form represents a combination of a bulk viscosity 
(linear in //y) and a von Neumann- Richtmyer viscosity (quadratic in /i^). 
The von Neumann-Richtmyer viscosity was initially introduced to suppress 
particle interpenetration in the presence of strong shocks. Eq. (17) provides a 
good treatment of shocks when a ~ 1, /3 ~ 2 and rf ~ 10 -2 (Monaghan 1989; 
Hernquist & Katz 1989). 



A well known problem with the classical AV of eq. (17) is that it can generate 
large amounts of spurious shear viscosity. For this reason, Hernquist & Katz 
(1989) introduced another form for the AV: 



n 



% + % if 



VI 



if (^ - vj) 



fj) < 
fj) > ' 



(19) 



where 

v\i + (3 Pl hj\V-v\^ if(V-tf). <0 




if (V • v), > (20) 
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and 



(V • v)i = — ^rrijivj - Vi) ■ V;Wy. 
Pi j 



(21) 



Although this form provides a slightly less accurate description of shocks than 
equation (17), it does exhibit less shear viscosity. 

More recently, Balsara (1995) has proposed the AV 




where 



(vi-Vj)-(ri-rj) fi+fj 



if (Vi - Vi) ■ (fi - fi) < 



if (vi - Vj) ■ {fi - fj) > 



Here fi is the form function for particle % defined by 



f. = (24) 

Jl \V-v\ i + \Vxv\ i + r ] 'c i /h i ' K } 



where the factor rj' ~ 10 4 — 10 5 prevents numerical divergences, (V • v)i is 
given by equation (21), and 

(V x v)i = - mj (vi -Vj)x V l W lJ . (25) 
Pi j 



The form function fi acts as a switch, approaching unity in regions of strong 
compression (|V • v\i » |V x v\i) and vanishing in regions of large vortic- 
ity (|V x v\i » |V • v\i). Consequently, this AV has the advantage that it 
is suppressed in shear layers. Note that since {pi/ p} +Pj/p]) ~ 2c 2 - / '(7 pif) , 
equation (22) behaves like equation (17) when |V • v\i » |V x v\i, provided 
one rescales the a and f3 in equation (22) to be a factor of 7/2 times the a 
and (3 in equation (17). 



1.2.4 Thermodynamics 

To complete the description of the fluid, either Ui or A4 is evolved according 
to a discretized version of the first law of thermodynamics. Although various 
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forms of these evolution equations exist, the most commonly used are 

^4z^-(;t+!+ n «) ® - ® ■ Vi ^-' (26) 



and 



dA 



7 



dt 2p? 



- H m i n *j - fy) ■ ViWij. 



(27) 



We call equation (26) the "energy equation," while equation (27) is the "en- 
tropy equation." Which equation one should integrate depends upon the prob- 
lem being treated. Each has its own advantages and disadvantages. Note that 
the derivation of equations (26) and (27) neglects terms proportional to the 
time derivative of hi. Therefore if we integrate the energy equation, even in the 
absence of AV, the total entropy of the system will not be strictly conserved if 
the particle smoothing lengths are allowed to vary in time; if the entropy equa- 
tion is used to evolve the system, the total entropy would then be strictly con- 
served when Ilij = 0, but not the total energy (Rasio 1991; Hernquist 1993). 
For more accurate treatments involving time- dependent smoothing lengths, 
see Nelson & Papaloizou (1994) and Serna et al. (1996). The energy equation 
has the advantage that other thermodynamic processes such as heating and 
cooling (Katz et al. 1996) and nuclear burning Garcia-Senz et al. 1998) can 
be incorporated more easily. 



1.2.5 Integration in Time 

The results of SPH simulations involving only hydrodynamic forces and gravity 
do not depend strongly on the actual time-stepping routine used, as long as the 
routine remains stable and accurate. A simple second-order explicit leap-frog 
scheme is often employed. Implicit schemes must be used when other processes 
such as heating and cooling are coupled to the dynamics (Katz et al. 1996). A 
low order scheme is appropriate for SPH because pressure gradient forces are 
subject to numerical noise. For stability, the timestep must satisfy a modified 
Courant condition, with hi replacing the usual grid separation. For accuracy, 
the timestep must be a small enough fraction of the dynamical time. 

Among the many possible choices for determining the timestep, the prescrip- 
tion proposed by Monaghan (1989) is recommended. This sets 

At = C N Mm(At u At 2 ), (28) 



where the constant dimensionless Courant number Cm typically satisfies 0.1 <> 
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Cat ^ 0.8, and where 



At x = Min, (hi/vi) 1 ' 2 , (29) 
At 2 = Mi^ ( — ^— — , ) , (30) 

with k being a constant of order unity. If the Hernquist & Katz AV [eq. (19)] 
is used, the quantity Maxj|/ijj| in equation (30) can be replaced by hi\V -v\i if 
(V • v)i < 0, and by otherwise. By accounting for AV-induced diffusion, the 
a and (3 terms in the denominator of equation (30) allow for a more efficient 
use of computational resources than simply using a smaller value of Cat. 



1.2.6 Smoothing Lengths and Accuracy 

The size of the smoothing lengths is often chosen such that particles roughly 
maintain some predetermined number of neighbors N^. Typical values of 
range from about 20 to 100. If a particle interacts with too few neighbors, 
then the forces on it are sporadic, a poor approximation to the forces on a 
true fluid element. In general, one finds that, for given physical conditions, the 
noise level in a calculation always decreases when Nn is increased. 

At the other extreme, large neighbor numbers degrade the resolution by re- 
quiring unreasonably large smoothing lengths. However, higher accuracy is 
obtained in SPH calculations only when both the number of particles iV and 
the number of neighbors N N are increased, with N increasing faster than 
N N so that the smoothing lengths hi decrease. Otherwise (e.g., if N is in- 
creased while maintaining constant) the SPH method is inconsistent, i.e., 
it converges to an unphysical limit. This can be shown easily by deriving the 
dispersion relation for sound waves propagating in simple SPH systems (Rasio 
1991). The choice of N N for a given calculation is therefore dictated by a com- 
promise between an acceptable level of numerical noise and the desired spatial 
resolution (which is ~ h oc l/N l J d in d dimensions) and level of accuracy. 



1.3 Results of Recent Test Calculations 



The authors and their collaborators have performed a series of systematic tests 
to evaluate the effects of spurious transport in SPH calculations. These tests 
are presented in detail in Lombardi et al. (1999), while here we summarize the 
main results. Our tests include (i) particle diffusion measurements, (ii) shock- 
tube tests, (iii) numerical viscosity measurements, and (iv) measurements of 
the spurious transport of angular momentum due to AV in differentially ro- 
tating, self-gravitating configurations. The results are useful for quantifying 
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the accuracy of the SPH scheme, especially for problems where shear flows 
or shocks are present, as well as for problems where true mixing is relevant. 
Other recent tests of SPH include those by Hernquist & Katz (1989) and by 
Steinmetz & Muller (1993). 



1.3.1 Particle Diffusion 

Many of our tests focus on spurious diffusion, the motion of SPH particles 
introduced as an artifact of the numerical scheme. Often applications require 
a careful tracing of the particle positions, and in these cases it is essential 
that spurious diffusion be small. For example, SPH simulations can be used 
to establish the degree of fluid mixing during stellar collisions, which is of 
primary importance in determining the subsequent stellar evolution of the 
merger remnants (see §2.1). It must be stressed that the amount of mixing 
determined by SPH calculations is always an upper limit. In particular, low- 
resolution calculations tend to be noisy, and this noise can lead to spurious 
diffusion of particles, independent of any real physical mixing of fluid elements. 

We have analyzed spurious diffusion by using SPH particles in a box with 
periodic boundary conditions to model a stationary fluid of infinite extent. 
For various noise levels (particle velocity dispersions) and neighbor numbers 
Nn, we measure the rate of diffusion, quantified by the diffusion coefficient 



Here the brackets () denote a time average, and Ar = (Arc 2 + Ay 2 + Az 2 ) 1 ^ 2 
is the total distance traveled by a particle due to spurious diffusion. Although 
strong shocks and AV in SPH calculations can lead to additional particle 
mixing (Monaghan 1989), particle diffusion is the dominant contribution to 
spurious mixing in weakly shocked fluids. 

Once expressed in terms of the number density of SPH particles and the sound 
speed, these diffusion coefficients can therefore be used to estimate spurious 
deviations in particle positions in a wide variety of applications, including 
self-gravitating systems. For each particle in some large-scale simulation, this 
spurious deviation is estimated simply by numerically integrating 

Ar 2 « [ Ddt. (32) 



The coefficient D in the integrand of equation (32) depends on the particle's 
velocity deviation from the local flow, the local number density n of particles, 
and the local sound speed c s , so that these quantities need to be monitored 
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for each particle during the simulation. Such a scheme was successfully used 
to estimate spurious mixing in the context of stellar collisions (Lombardi et 
al. 1996), where typically (with N = 3 x 10 4 and N N m 64) the diffusion 
coefficient was very roughly D ~ 0.05c s n -1 / 3 . 

For sufficiently low noise levels, the diffusion coefficient essentially vanishes, 
as the particles simply oscillate around equilibrium lattice sites. We say that 
such a system has "crystallized." For a neighbor number Nn ps 64, a system 
of SPH particles will crystallize if the root mean square velocity dispersion 
is less than about 3-4% of the sound speed. We find that crystallized cu- 
bic lattices are unstable against perturbations, while lattice types with large 
packing fractions, such as hexagonal close-packed, are stable. For this reason 
it may sometimes be better to construct initial data by placing particles in an 
hexagonal close-packed lattice, rather than in a cubic lattice as is often done. 

The diffusion coefficients have been measured using equal-mass particles. Some- 
times, however, SPH simulations use particles of unequal mass so that less 
dense regions can still be highly resolved. To test the effects of unequal mass 
particles in a self-gravitating system, we constructed an equilibrium n — 1.5 
polytrope (a polytrope is an idealized model for a spherical star, characterized 
by a relation of the form P = p 1 between pressure P and density p; the poly- 
tropic index n is defined by 7 = 1 + 1/n), using particle masses which increased 
with radius in the initial configuration. Allowing the system to evolve, we ob- 
served that the heaviest particles gradually migrated towards the center of the 
star, exchanging places with less massive particles. For a polytrope modeled 
with N 1.4 x 10 4 particles and a neighbor number Nn ~ 64, the distribution 
of particle masses is reversed within roughly 80 dynamical timescales. This is 
caused by the interactions among neighboring particles via the smoothing 
kernel. These interactions allow energy exchange, and equipartition of energy 
then requires the heavier particles to sink into the gravitational potential well. 
Spurious mixing is therefore a more complicated process in simulations which 
use unequal mass particles: each particle has a preferred direction to migrate, 
and in a dynamical application this direction can be continually changing. For 
simulations in which fluid mixing is important, equal- mass particles are an 
appropriate choice. 



1.3.2 Shock Tube Tests 

The diffusion tests just described are all done in the absence of shocks and 
without AV. To test the AV schemes described in §1.2, we turn to a periodic 
version of the 1-D Riemann shock-tube problem. Initially, fluid slabs with con- 
stant (and alternating) density p and pressure p are separated by an infinite 
number of planar, parallel, and equally spaced interfaces. We treat this in- 
herently 1-D problem with both a 1-D and a 3-D SPH code. The 1-D code is 
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naturally more accurate, and provides a benchmark against which we can com- 
pare the results of our 3-D code. In both cases, periodic boundary conditions 
allow us to model the infinite number of slabs. 

Using various values of a and (3, we performed a number of such shock tube cal- 
culations with our 3-D code, at both Mach numbers Ai ~ 1.6 and M. 13.2. 
We then compared the time variation of the internal energy and entropy of 
the system against that of the 1-D simulation. Furthermore, since any motion 
perpendicular to the bulk fluid flow is spurious, we were also able to examine 
spurious mixing in these simulations. We find that all three forms of AV can 
handle shocks well. For example, with N = 10 4 and Nn ~ 64, there is better 
than 2% agreement with the 1-D code's internal energy vs. time curve when 
M. ~ 1.6, and agreement at about the 3% level when M. m 13.2. We also 
find that both equations (17) and (22), as compared to equation (19), allow 
less spurious mixing and do somewhat better at reproducing the 1-D code's 
results. 

Such simulations are a useful and realistic way to calibrate spurious transport, 
since the test problem, which includes shocks and significant fluid motion, has 
many of the same properties as real astrophysical problems. In fact, the recoil 
shocks in stellar collisions do tend to be nearly planar, so that even the 1-D 
geometry of the shock fronts is realistic. The periodic boundary conditions 
play the role of gravity in the sense that they prevent the gas from expanding 
to infinity. 



1.3.3 Shear Flows 

To test the various AV forms in the presence of a shear flow, we impose the 
so-called slipping boundary conditions on a periodic box, as is commonly done 
in molecular dynamics (see, e.g., Naitoh & Ono 1976). The resulting "station- 
ary Couette flow" has a velocity field close to (v x ,v y ,v z ) = (vq]j/L, 0,0) and 
allows us to measure the numerical viscosity of the particles. As in the shock 
tube tests, we also examine spurious mixing in the direction perpendicular to 
the fluid flow. These shear tests therefore allow us to further investigate the 
accuracy of our SPH code as a function of the AV parameters and scheme. 
We find that both the Hernquist & Katz AV [eq. (19)] and the Balsara AV 
[eq. (22)] exhibit less viscosity than the classical AV [eq. (17)]. However, the 
classical AV does allow significantly less spurious mixing than the other forms. 
For all three forms of the AV, increasing a and (3 tends to damp out the noise 
and consequently decrease spurious mixing, but it also increases the spurious 
shear viscosity. 

Rotation plays an important role in many hydrodynamic processes. For in- 
stance, a collision between stars can yield a rapidly and differentially rotating 
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merger remnant. Even in the absence of shocks, AV tends to damp away differ- 
ential rotation due to the relative velocity of neighboring particles at slightly 
different radii, and an initially differentially rotating system will tend towards 
rigid rotation on the viscous dissipation timescale. In systems best modeled 
with a perfect fluid, ideally with a viscous timescale r = oo, any such angular 
momentum transport introduced by the SPH scheme is spurious. 

As a concrete example, we consider an axisymmetric equilibrium configuration 
differentially rotating with an angular velocity profile Vi{w) oc -oj~ x , where w is 
the distance from the rotation axis and A is a constant of order unity. We then 
analytically estimate the viscous dissipation timescale for each of the three 
AVs discussed in §1.2. These analytic estimates are found to closely match 
numerically measured values of the timescale. Both the Hernquist & Katz AV 
[eq. (19)] and the Balsara AV [eq. (22)] yield longer viscous timescales than 
the classical AV [eq. (17)], and hence are better at maintaining the angular 
velocity profile. The Balsara AV clearly does best in this regard, with a viscous 
timescale roughly N]J 2 times larger than for the classical AV. 

When choosing values of AV parameters, one must weigh the relative impor- 
tance of shocks, shear, and fluid mixing. For this reason, it is an application- 
dependent, somewhat subjective matter to specify "optimal values" of a and 
(3. We do, however, roughly delineate the boundaries of the region in parameter 
space that gives acceptable results in Lombardi et al. (1999). 

Our results concerning the various AV forms can be summarized as follows 
(see Lombardi et al. 1999 for more details). We find that the AVs defined by 
equations (17) and (22) do equally well both in their handling of shocks and 
in their controlling of spurious mixing, and do slightly better than equation 
(19). Furthermore, both equations (19) and (22) do introduce less numerical 
viscosity than equation (17). Since equation (22), Balsara's form of AV, does 
indeed significantly decrease the amount of shear viscosity without sacrificing 
accuracy in the treatment of shocks, we conclude that it is an appropriate 
choice for a broad range of problems. This is consistent with the successful 
use of Balsara's AV reported by Navarro & Steinmetz (1997) in their models 
of rotating galaxies. 



2 SPH Calculations of Stellar Interactions 

The vast majority of recent 3-D calculations of dynamical interactions between 
stars have been done using the SPH method. These include collisions (Benz & 
Hills 1992; Lombardi et al. 1996), binary coalescence (Davies et al. 1994; Rasio 
& Shapiro 1995), common envelope evolution (Rasio & Livio 1996; Terman & 
Taam 1996), accretion flows (Bate & Bonnell 1997; Theuns et al. 1996), and 
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tidal disruption (Laguna et al. 1993). 

SPH has many advantages over more traditional methods of numerical hy- 
drodynamics for these calculations. First, and perhaps most importantly, the 
advection of the fluid while the stars are moving along their initial trajectories 
is handled very easily by SPH. For example, in the case of binary coalescence 
(see §2.2 below), one often has to follow the motion of the two stars for several 
orbital periods before the final merger occurs. Merely tracking the motion of a 
star across a large 3-D grid for many dynamical times can be very challenging 
when using an Eulerian scheme, especially in the presence of a sharp stellar 
surface (as in the case of neutron stars, which contain a fairly incompressible 
fluid). In addition, other physical processes can be studied much more easily 
using a Lagrangian scheme. An example is hydrodynamic mixing, which is 
a crucial process in the study of certain stellar merger processes (see §2.1). 
Since chemical abundances are passively advected quantities during a dynam- 
ical evolution, the chemical composition in the final fluid configuration can be 
determined after the completion of a calculation simply by noting the original 
and final positions of all SPH particles and by assigning particle abundances 
according to an initial profile. The adaptiveness of the scheme, with particles 
automatically concentrating in regions of higher density is also an important 
advantage, although this is more crucial in situations involving large density 
contrasts, such as in simulations of star formation or galaxy formation. 

2.1 Stellar Collisions 

As a first illustration of the use of SPH for a typical stellar interaction problem, 
we summarize in this section recent work on the numerical calculation of 
collisions between two main-sequence stars. 

2.1.1 Motivation 

Close dissipative encounters and direct physical collisions between stars occur 
frequently in dense star clusters. The dissipation of kinetic energy in close stel- 
lar encounters can have a direct influence on the overall dynamical evolution 
of a cluster. Observational evidence for stellar collisions and mergers in glob- 
ular clusters is provided by the existence of large numbers of blue stragglers 
in these systems. These are peculiar main-sequence (hydrogen burning) stars 
that appear younger and more massive than all other, normal main-sequence 
stars in the cluster. 

Blue stragglers have long been thought to be formed through the merger of 
two lower-mass stars, either in a collision or following binary coalescence (see, 
e.g., the review by Livio 1993). Clear indication for a collisional origin of blue 
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Fig. 1. Snapshots of density contours in the orbital plane for a parabolic collision 
between two main-sequence stars of masses Mi and M2 = 0.75M\. The impact 
parameter has been chosen such that the corresponding point-mass orbit would 
have a pericenter separation r p = 0.25(i?i + R2), where R\ and R2 = 0.56i?i are 
the stellar radii. There are eight density contours, which are spaced logarithmically 
and cover four decades down from the maximum. The elapsed time in the upper left 
corner of each frame is in units of the dynamical timescale (Rf/GMi) 1 / 2 . Adapted 
from Lombardi et al. (1996). 

stragglers has come from observations of globular clusters with the Hubble 
Space Telescope. Large numbers of blue stragglers were found to be concen- 
trated in the cores of the densest clusters (such as M15 and M30). 

Following early numerical work in 2-D (e.g., Shara & Shaviv 1978), Benz 
& Hills (1987, 1992) performed the first 3-D calculations of direct collisions 
between two main sequence stars using SPH. An important result of this 
pioneering study was that collisions could lead to a thoroughly mixed merger 
remnant. The mixing of fresh hydrogen fuel into the core of the remnant could 
then reset its nuclear clock, allowing the blue straggler to burn hydrogen for 
a full main-sequence lifetime (tMS ~ 10 9 yr) after its formation. 



2.1.2 Recent Results using SPH 

The authors and their collaborators have re-examined collisions between main- 
sequence stars, and, in particular, the question of mixing during mergers, by 
performing a set of numerical hydrodynamic calculations using SPH (Lom- 
bardi, Rasio, & Shapiro 1995, 1996; Sills et al. 1997; Sills & Lombardi 1997). 
This new work differs from the previous study of Benz & Hills (1987) by 
adopting more realistic models of globular cluster stars, and by performing 
numerical calculations with increased spatial resolution. 

Snapshots from one of our recent calculations are shown in Figure 1. The 
main new results of these SPH calculations can be summarized as follows. 
Merger remnants produced by parabolic collisions are always far from chemi- 
cally homogeneous. In the case of collisions between two nearly identical stars, 
the amount of hydrodynamic mixing during the collision is minimal. In fact, 
the final chemical composition profile is very close to the initial profile of the 
parent stars. For two turnoff stars (i.e., close to hydrogen exhaustion at the 
center), this means that the merger remnant is born with very little hydro- 
gen to burn in its core and, consequently, that the object may not be able to 
remain on the main sequence for long. In the case of a collision between two 
stars of different masses, the chemical composition profile of the merger rem- 
nant tends to be more homogeneous, but it remains true that little hydrogen 
is injected into its core. 
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At a qualitative level, these results can be understood very simply in terms 
of the requirement of convective (dynamical) stability of the final hydrostatic 
equilibrium configurations. For non-rotating remnants, convective stability re- 
quires that the entropic variable A [see eq. (1)] increase monotonically from 
the center to the surface (the so-called Ledoux criterion). If shock- heating 
could be neglected entirely (which is not unreasonable for the low-velocity 
collisions occurring in globular clusters), then one could predict the final com- 
position profile of a merger remnant simply by observing the composition and 
A profiles of the parent stars. Since fluid elements conserve both their value 
of A (in the absence of shocks) and their composition, the final composition 
profile of a remnant is constructed simply by combining mass shells in order 
of increasing A, from the center to the outside. Many features of the results 
follow directly. For example, in the case of a collision between two identical 
stars, it is obvious why the composition profile of the merger remnant remains 
very similar to that of the parent stars, since shock-heating is significant only 
in the outer layers of the stars, which contain a very small fraction of the total 
mass. Furthermore, the dense, helium-rich material is concentrated in the deep 
interior of the parent stars, where shock-heating is negligible, and therefore it 
remains concentrated in the deep interior of the final configuration. 

2.2 Coalescing Compact Binaries 

As a second illustration of the use of SPH for calculations of stellar interac- 
tions, we now turn to the dynamical evolution of compact binary star sys- 
tems (containing two compact objects — black holes or neutron stars — in orbit 
around one another). 

2.2.1 Motivation 

Coalescing compact binaries are the most promising known sources of gravi- 
tational radiation that could be detected by the new generation of laser in- 
terferometers now under construction. These include the Caltech-MIT LIGO 
(Abramovici et al. 1992; Cutler et al. 1993) and the European projects VIRGO 
(Bradaschia et al. 1990) and GEO 600 (Danzmann 1998). In addition to 
providing a major new confirmation of Einstein's theory of general relativ- 
ity, including the first direct proof of the existence of black holes (Flanagan 
& Hughes 1998; Lipunov et al. 1997), the detection of gravitational waves 
from coalescing binaries at cosmological distances could provide accurate in- 
dependent measurements of the Hubble constant and mean density of the 
universe (Schutz 1986; Chernoff & Finn 1993; Markovic 1993). For an ex- 
cellent recent review on the detection and sources of gravitational radiation, 
see Thorne (1996). Coalescing compact binaries are important for other areas 
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of astrophysics as well. In particular, many theoretical models of gamma-ray 
bursts have postulated that the energy source for the bursts could be coalesc- 
ing compact binaries at cosmological distances (Eichler et al. 1989; Narayan, 
Paczyhski, & Piran 1992). For a discussion of the hydrodynamics of binary 
coalescence in the context of gamma-ray burst models, see Ruffert et al. (1997) 
and references therein. 

Here we focus on binaries containing two neutron stars (hereafter NS), for 
which the final coalescence is a purely hydrodynamic process that has been 
well-studied using SPH (as well as other methods). 



2.2.2 Hydrodynamics of the Binary Merger Process 

Hydrostatic equilibrium configurations for binary systems with sufficiently 
close components can become dynamically unstable (Chandrasekhar 1975; Tas- 
soul 1975). The physical nature of this instability is common to all binary in- 
teraction potentials that are sufficiently steeper than 1/r (see, e.g., Goldstein 
1980, §3.6). It is analogous to the familiar instability of test particles in circu- 
lar orbits sufficiently close to a black hole (Shapiro & Teukolsky 1983, §12.4). 
Here, however, it is the tidal interaction that is responsible for the steepening 
of the effective interaction potential between the two stars and for the desta- 
bilization of the circular orbit. The physical properties of this instability, and 
its consequences for NS binary coalescence, have been studied by the authors 
and their collaborators, both numerically using SPH (Rasio & Shapiro 1994, 
1995, hereafter RS; Rasio 1998) and using analytical perturbation methods 
(Lai et al. 1993, 1994, hereafter LRS; Lombardi et al. 1997). 

The instability was first identified by RS using SPH calculations where the 
evolution of binary equilibrium configurations was calculated for two identical 
polytropes with 7 = 2. It was found that when r <, 3R (r is the binary sepa- 
ration and R the radius of an unperturbed NS), the orbit becomes unstable to 
radial perturbations and the two stars undergo rapid coalescence. For r ^> 3R, 
the system could be evolved dynamically for many orbital periods without 
showing any sign of orbital evolution (in the absence of dissipation). Many 
of the results derived in RS and LRS concerning the stability properties of 
NS binaries have been confirmed recently in completely independent work by 
New & Tohline (1997) and by Zhuge, Centrella, & McMillan (1996). New & 
Tohline (1997) used completely different numerical methods (a combination 
of a 3-D Self-Consistent Field code for constructing equilibrium configurations 
and a grid-based Eulerian code for following the dynamical evolution of the 
binaries), while Zhuge et al. (1996) used SPH. 
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Fig. 2. Evolution of an unstable binary containing two identical neutron stars. The 
stars are modeled as polytropes with 7 = 3, corresponding to a stiff nuclear equation 
of state. Projections of all SPH particles onto the orbital plane of the binary are 
shown at different times (the orbital motion is in the counterclockwise direction). 
Units are such that G = M = R = 1, where G is the gravitational constant and 
R and M are the radius and mass of an unperturbed (spherical) star. The initial 
orbital period is P or & ~ 24 in these units. Adapted from Rasio & Shapiro (1994). 



2. 2. 3 Typical SPH Results 

For simplicity, we describe here the dynamical evolution of an unstable, ini- 
tially synchronized binary containing two identical stars. Typical SPH results 
for this case are shown in Figure 2. During the initial, linear stage of the insta- 
bility, the two stars approach each other and come into contact after about one 
orbital revolution. In the corotating frame of the binary, the relative velocity 
remains very subsonic, so that the evolution is adiabatic at this stage. This 
is in sharp contrast to the case of a head-on collision between two stars on a 
free-fall, radial orbit, where shocks can be very important for the dynamics. 
Here the stars are constantly being held back by a (slowly receding) centrifugal 
barrier, and the merging, although dynamical, is much more gentle. 

After typically 2-3 orbital revolutions the innermost cores of the two stars 
have merged and the system resembles a single, very elongated ellipsoid. At 
this point a secondary instability occurs: mass shedding sets in rather abruptly. 
Material is ejected through the outer Lagrange points of the effective potential 
and spirals out rapidly. In the final stage, the spiral arms widen and merge 
together. The relative radial velocities of neighboring arms as they merge are 
supersonic, leading to some shock-heating and dissipation. As a result, a hot, 
nearly axisymmetric rotating halo forms around the central dense core. The 
halo contains about 20% of the total mass and the rotation profile is close to 
a pseudo-barotrope (Tassoul 1978), with the angular velocity decreasing as a 
power-law fl oc w~ v where u <, 2 and w is the distance to the rotation axis. 
The core is rotating uniformly near breakup speed and contains about 80% 
of the mass still in a cold, degenerate state. If the initial NS had masses close 
to 1.4 Mq, then most recent stiff EOS would predict that the final merged 
configuration is still stable and will not immediately collapse to a black hole, 
although it might ultimately collapse to a black hole as it continues to lose 
angular momentum. The emission of gravitational radiation during dynamical 
coalescence can be calculated perturbatively using the quadrupole approxima- 
tion (RS). Both the frequency and amplitude of the emission peak somewhere 
during the final dynamical coalescence, typically just before the onset of mass 
shedding. Immediately after the peak, the amplitude drops abruptly as the 
system evolves towards a more axially symmetric state. 



18 



References 

[I] Abramovici, M., et al. 1992, Science, 256, 325 
[2] Balsara, D. 1995, J. Comp. Phys., 121, 357 

[3] Bate, M.R., & Bonnell, I. A. 1997, MNRAS, 285, 33 

[4] Benz, W., & Hills, J. G. 1987, ApJ, 323, 614 

[5] Benz, W., & Hills, J. G. 1992, ApJ, 389, 546 

[6] Bradaschia, C, et al. 1990, Nucl. Instr. Methods A, 289, 518 

[7] Burkert, A., Bate, M.R., & Bodenheimer, P. 1997, MNRAS, 289, 497 

[8] Chandrasekhar, S. 1975, ApJ, 202, 809 

[9] Chernoff, D.F., & Finn, L.S. 1993, ApJ, 411, L5 

[10] Couchman, H.M.P, Thomas, PA., & Pearce, F.R. 1995, ApJ, 452, 797 

[II] Cutler, C, et al. 1993, PRL, 70, 2984 

[12] Danzmann, K. 1998, in Relativistic Astrophysics, eds. H. Riffert et al. (Proc. 
of 162nd W.E. Heraeus Seminar, Wiesbaden: Vieweg Verlag), 48 

[13] Dave, R., Dubinski, J., & Hernquist, L. 1997, New A, 2, 277 

[14] Davies, M.B., Benz, W., Piran, T., & Thielemann, F.K. 1994, ApJ, 431, 742 

[15] Eichler, D., Livio, M., Piran, T., & Schramm, D.N. 1989, Nature, 340, 126 

[16] Evrard, A.E. 1988, MNRAS, 235, 911 

[17] Flanagan, E.E., & Hughes, S.A. 1998, PRD, 57, 4566 

[18] Garcia-Senz, D., Bravo, E., & Serichol, N. 1998, ApJS, 115, 119 

[19] Gingold, R.A., & Monaghan, J.J. 1977, ApJ, 181, 375 

[20] Gingold, R.A., & Monaghan, J.J. 1982, J. Comp. Phys., 46, 429 

[21] Goldstein, H. 1980, Classical Mechanics (Reading: Addison-Wesley) 

[22] Herant, M., & Benz, W. 1992, ApJ, 387, 294 

[23] Hernquist, L. 1993, ApJ, 404, 717 

[24] Hernquist, L., & Katz, N. 1989, ApJS, 70, 419 

[25] Katz, N. 1992, ApJ, 391, 502 

[26] Katz, N., Weinberg, D.H., & Hernquist, L. 1996, ApJS, 105, 19 
[27] Klessen, R. 1997, MNRAS, 292, 11 



19 



[28] Laguna, P., Miller, W.A., Zurek, W.H., k Davies, M.B. 1993, ApJ, 410, L83 
[29] Lai, D., Rasio, F.A., k Shapiro, S.L. 1993, ApJS, 88, 205 
[30] Lai, D., Rasio, F.A., k Shapiro, S.L. 1994, ApJ, 437, 742 

[31] Livio, M. 1993, in ASP Conf. Ser. Vol. 53, Blue Stragglers, ed. R. A. Saffer (San 
Francisco: ASP), 3 

[32] Lombardi, J.C., Jr., Rasio, F.A., k Shapiro, S.L. 1995, ApJ, 445, L117 

[33] Lombardi, J.C., Jr., Rasio, F.A., k Shapiro, S.L. 1996, ApJ, 468, 797 

[34] Lombardi, J.C., Jr., Rasio, F.A., k Shapiro, S.L. 1997, PRD, 56, 3416 

[35] Lombardi, J.C., Jr., Rasio, F.A., Sills, A. k Shapiro, S.L. 1999, J Comp Phys, 
in press 

[36] Lucy, L. B. 1977, AJ, 82, 1013 

[37] Markovic, D. 1993, PRD, 48, 4738 

[38] Monaghan, J.J. 1985, Comp. Phys. Rep., 3, 71 

[39] Monaghan, J.J. 1989, J. Comp. Phys., 82, 1 

[40] Monaghan, J.J. 1992, ARA&A, 30, 543 

[41] Monaghan, J.J., k Lattanzio, J.C. 1985, A&A, 149, 135 

[42] Naitoh, T., k Ono, S. 1976, Physics Letters, 57A, 448 

[43] Narayan, R., Paczyhski, B., k Piran, T. 1992, ApJ, 395, L83 

[44] Navarro, J.F., k Steinmetz, M. 1997, 478, 13 

[45] Nelson, A.F., Benz, W., Adams F.C., k Arnett, D. 1998, ApJ, in press 

[46] Nelson, R.P., k Papaloizou, J.C.B. 1994, MNRAS, 270, 1 

[47] New, K.C.B., k Tohline, J.E. 1997, ApJ, 490, 311 

[48] Phinney, E.S. 1991, ApJ, 380, L17 

[49] Rasio, F.A. 1991, PhD Thesis (Cornell University) 

[50] Rasio, F.A. 1995, ApJ, 444, L41 

[51] Rasio, F.A. 1998, in The Many Faces of Neutron Stars, eds. J. van Paradijs et 
al. (Dordrecht: Kluwer), in press 

[52] Rasio, F.A., k Livio, M. 1996, ApJ, 471, 366 

[53] Rasio, F.A., k Shapiro, S.L. 1992, ApJ, 401, 226 

[54] Rasio, F.A., k Shapiro, S.L. 1994, ApJ, 432, 242 

[55] Rasio, F.A., k Shapiro, S.L. 1995, ApJ, 438, 887 



20 



[56] Ruffert, M., Janka, H.-T., Takahashi, K., & Schafer, G. 1997, A&A, 319, 122. 
[57] Schutz, B.F. 1986, Nature, 323, 310. 

[58] Serna, A., Alimi, J.-M., & Chieze, J.-P. 1996, ApJ, 461, 884 

[59] Shapiro, S.L., & Teukolsky, S.A. 1983, Black Holes, White Dwarfs, and Neutron 
Stars (New York: Wiley) 

[60] Shapiro, P.R., Martel, H., Villumsen, J.V., Owen, J.M. 1996, ApJS, 103, 269 

[61] Shara, M.M., & Shaviv, G. 1978, MNRAS, 183, 687 

[62] Sills, A., & Lombardi, J.C., Jr. 1997, ApJ, 484, L51 

[63] Sills, A., Lombardi, J.C., Jr., Bailyn, C.D., Demarque, P., Rasio, F.A., & 
Shapiro, S.L. 1997, ApJ, 487, 290 

[64] Steinmetz, M. 1996, MNRAS, 278, 1005 

[65] Steinmetz, M., & Miiller, E. 1993, A&A, 268, 391 

[66] Tassoul, M. 1975, ApJ, 202, 803 

[67] Tassoul, J. 1978, Theory of Rotating Stars (Princeton: Princeton Univ. Press) 

[68] Terman, J.L., & Taam, R.E. 1996, ApJ, 458, 692 

[69] Theuns, T., Boffin, H.M.J. , & Jorissen, A. 1996, MNRAS, 280, 1264 

[70] Thorne, K.S. 1996, in Compact Stars in Binaries, IAU Symp. 165, eds. J. van 
Paradijs et al. (Kluwer, Dordrecht), 153 

[71] Will, C.M. 1994, in Relativistic Cosmology, ed. M. Sasaki (Universal Academy 
Press), 83 

[72] Zhuge, X., Centrella, J.M., & McMillan, S.L.W. 1996, PRD, 54, 7261 



21 



